function [ Supp ] = RunOMP_forMB_SPECIAL(Y,A,N, MoreIters,y)
[mA,nA] = size(A);
[mY,nY] = size(Y);

Fnorm=[];
Fnorm2=[];
residual = Y;
OneSidedSupp = [];
Supp = [];
SymmetricSupp = [];
iter = 1;

NormACols = sqrt(diag(A'*A));

    while (iter <= N+MoreIters)
        if  (iter <= N)
            % Matching step
            Z_1 = A'*residual;
            Z = sqrt(sum(abs(Z_1).^2,2))./NormACols;
            [maxVal, maxPos] = max(Z);
            BestLoc = maxPos(1);
            % Update support
            OneSidedSupp = [OneSidedSupp BestLoc];
            SymmetricLoc =  (nA+1-BestLoc);
            if (BestLoc ~= SymmetricLoc)
                SymmetricSupp = [SymmetricSupp SymmetricLoc];
            end
            Supp = [OneSidedSupp  SymmetricSupp];
        end
    %     if (iter == N)          %  Now we search for the adjacent positions
    %         CurrentLen = length(OneSidedSupp);
    %         for i=1:CurrentLen
    %             DominantBand = OneSidedSupp(i);
    %             CheckSupp = [DominantBand-1   DominantBand+1];
    %             CheckSupp = min(CheckSupp,nA);
    %             CheckSupp = max(CheckSupp,1);
    %             CheckSupp = unique(CheckSupp);
    %             if  (   (not(isempty(intersect(CheckSupp,OneSidedSupp))))   || ...
    %                     (not(isempty(intersect(CheckSupp,SymmetricSupp)))) )
    %                 % meaning that we already picked the adjacent position
    %                 continue;
    %             end
    %             % Matching step
    %             Z_2 = A(:,CheckSupp)'*residual;
    %             Z = sqrt(sum(abs(Z_2).^2,2))./NormACols(CheckSupp);
    %             [maxVal, maxPos] = max(Z);
    %             BestLoc = CheckSupp(maxPos(1));
    %             % Update support
    %             OneSidedSupp = [OneSidedSupp BestLoc];
    %             SymmetricLoc =  (nA+1-BestLoc);
    %             if (BestLoc ~= SymmetricLoc)
    %                 SymmetricSupp = [SymmetricSupp SymmetricLoc];
    %             end
    %             Supp = [OneSidedSupp  SymmetricSupp];
    %         end
    %     end
        % Project residual
        As = A(:,Supp);
        solution = As*pinv(As)*Y;
        residual = Y-solution;
        R=(residual')*residual;
%         Fnorm(iter) = sum(sqrt(diag(R))) / nY ;
        
%         %% debug
%         solution2 = As*pinv(As)*y;
%         residual2= y-solution2;
%         R2=(residual2')*residual2;
%         Fnorm2(iter) = sum(sqrt(diag(R2))) / length(y) ;
        
        %%
        iter = iter+1;

    end
% Fnorm
% Fnorm2
%% DEBUG
% R=(residual')*residual;
% Fnorm = sum(diag(R));
% g=6;

